      Program dGTSTHarmonic
C***********************************************************************
C     This program solves the dGT equation including ST
C     The parameters are set for the Lar expt.
C***********************************************************************
      Parameter(n=1024)
      Complex*16 A,ii,part
      Real*8 h,pi,s,exptdata,delta,epsilon,k0,omega0,gamma
      Integer j,q,carloc
      Dimension A(0:n-1),exptdata(1:3,1:8192)
      Common A
      pi=4.0d0*datan(1.0d0)
      ii=dcmplx(0.0d0,1.0d0)

      open(unit=1,file='ICST.out',status='old')
      read(1,*) exptdata
      close(unit=1)

      carloc=157
      omega0=exptdata(1,carloc)
      k0=1.52759556257556977d0
      epsilon=2.d0*k0*dsqrt(exptdata(2,carloc)**2+exptdata(3,carloc)**2)
      delta=7.229992843192701d0
      s=omega0*2.0d0*pi*epsilon*exptdata(1,2)**(-1)/2.0d0
      gamma=k0**2*71.70d0/981.0d0
      h=k0*epsilon**2/10.0d0
      q=5000

      Do 12 j=0,n-1
         A(j)=0.0d0
 12   Continue

      Do 15 j=0,n-1
         Do 11 jj=-39,39
            part=exptdata(2,carloc+jj)+ii*exptdata(3,carloc+jj)
            A(j)=A(j)+part*exp(jj*ii*2.0d0*pi*j/(n*1.0d0))
 11      Continue
 15   Continue

      Do 13 j=0,n-1
         A(j)=k0/epsilon*A(j)
 13   Continue

      open(unit=2,file="dGTSTHarmonicAlldata.out")
      open(unit=3,file="dGTSTHarmonicAllCQ.out")

      s=s*1.0d0/pi
      Call sixthorder(q,s,h,k0,epsilon,delta,gamma)

      close(2)
      close(3)

      End



      Subroutine linear(s,timestep,k0,delta,epsilon,gamma)
C***********************************************************************
C     This subroutine solves the linear part of vDysthe.
C***********************************************************************

      Complex*16 A,ii
      Parameter(n=1024)
      Integer i
      Real*8 l,s,timestep,k0,AR,delta,epsilon,gamma
      Real*8 mu1,mu3,mu7,mu8
      Dimension A(0:n-1),AR(0:1,0:n-1)
      Common A
      Equivalence(A,AR)
      ii=dcmplx(0.0d0,1.0d0)

      Call FFT(AR(0,0),AR(1,0),2,n,1)

      mu1=(1.0d0+gamma)*(1.0d0-6.0d0*gamma-3.0d0*gamma**2)
      mu1=mu1/(1.0d0+3.0d0*gamma)**3
      mu3=4.0d0*gamma*(1.0d0+gamma)**2
      mu3=mu3*(5.0d0-10.0d0*gamma-3.0d0*gamma**2)
      mu3=mu3/(1.0d0+3.0d0*gamma)**5
      mu7=(5.0d0+10.0d0*gamma+9.0d0*gamma**2)/(1.0d0+3.0d0*gamma)**2
      mu8=2.0d0*(1.0d0+gamma)*(-5.00+10.0d0*gamma+3.0d0*gamma**2)
      mu8=mu8/(1.0d0+3.0d0*gamma)**3

      Do 10 i=0,n/2-1
         l=-1.0d0*i/s
         A(i)=A(i)*exp(-mu1*ii*l*l*2.0d0*timestep)
         A(i)=A(i)*exp(epsilon*mu3*ii*l*l*l*2.0d0*timestep)
         A(i)=A(i)*exp(-delta*2.0d0*timestep)
         A(i)=A(i)*exp(-mu7*epsilon*delta*l*2.0d0*timestep)
         A(i)=A(i)*exp(mu8*epsilon**2*delta*l*l*2.0d0*timestep)
 10   Continue
      
      Do 20 i=n/2+1,n-1
         l=1.0d0*(n-i)/s
         A(i)=A(i)*exp(-mu1*ii*l*l*2.0d0*timestep)
         A(i)=A(i)*exp(epsilon*mu3*ii*l*l*l*2.0d0*timestep)
         A(i)=A(i)*exp(-delta*2.0d0*timestep)
         A(i)=A(i)*exp(-mu7*epsilon*delta*l*2.0d0*timestep)
         A(i)=A(i)*exp(mu8*epsilon**2*delta*l*l*2.0d0*timestep)
 20   Continue
      




C**** Notice the 3/8ths rule!!!
      Do 22 i=3*n/8,5*n/8
         A(i)=0.0d0
 22   Continue

      Call FFT(AR(0,0),AR(1,0),2,n,-1)
      
      Do 30 i=0,n-1
         A(i)=A(i)/(1.0d0*n)
 30   Continue

      Return
      End

      Subroutine nonlinear(timestep,k0,epsilon,s,gamma)
C***********************************************************************
C       This subroutine solves the nonlinear part of vDysthe.
C***********************************************************************

      Complex*16 A,ii,k1,k2,k3,k4
      Real*8 timestep,deltattt,k0,epsilon,s,gamma
      Integer i
      Parameter(n=1024)
      Dimension A(0:n-1),k1(0:n-1),k2(0:n-1),k3(0:n-1),k4(0:n-1)
      Common A
      ii=dcmplx(0.0d0,1.0d0)

      deltattt=timestep/1.0d0

      Do 20 i=1,1

         Call RHS(A,k0,epsilon,k1,s,gamma)
         Call RHS(A+deltattt*k1,k0,epsilon,k2,s,gamma)
         Call RHS(A+deltattt*k2,k0,epsilon,k3,s,gamma)
         Call RHS(A+2.0d0*deltattt*k3,k0,epsilon,k4,s,gamma)
         A=A+deltattt/3.0d0*(k1+2*(k2+k3)+k4)

 20   Continue

      End


      Subroutine RHS(B,k0,epsilon,k,s,gamma)
C***********************************************************************
C***********************************************************************

      Complex*16 A,B,ii,k,Bx,Phix
      Real*8 k0,epsilon,BxR,s,PhixR,gamma,mu2,mu4,mu5,mu6
      Integer i,x
      Parameter(n=1024)
      Dimension A(0:n-1),B(0:n-1),k(0:n-1),Bx(0:n-1),BxR(0:1,0:n-1)
      Dimension Phix(0:n-1),PhixR(0:1,0:n-1)
      Common A
      Equivalence(Bx,BxR)
      Equivalence(Phix,PhixR)
      ii=dcmplx(0.0d0,1.0d0)

      mu2=(8.0d0+gamma+2.0d0*gamma**2)
      mu2=mu2/(2.0d0+2.0d0*gamma-12.0d0*gamma**2)
      mu4=(1.0d0+gamma)*(1.0d0-3.0d0*gamma)*(8.0d0+gamma+2.0d0*gamma**2)
      mu4=mu4/(-1.0d0+2.0d0*gamma)/(1.0d0+3.0d0*gamma)**3
      mu5=2.0d0*(1.0d0+gamma)
      mu5=mu5*(-16.0d0+13.0d0*gamma-55.0d0*gamma**2+12.0d0*gamma**4)
      mu5=mu5/(1.0d0-2.0d0*gamma)**2/(1.0d0+3.0d0*gamma)**3
      mu6=-16.0d0*(1.0d0+gamma)**3/(1.0d0+3.0d0*gamma)**3

      Do 10 x=0,n-1
         Bx(x)=B(x)
 10   Continue

      Call FFT(BxR(0,0),BxR(1,0),2,n,1)
      Do 20 i=0,n/2
         Bx(i)=-ii*i/s*Bx(i)
 20   Continue
      Do 30 i=n/2+1,n-1
         Bx(i)=ii*(n-i)/s*Bx(i)
 30   Continue
      Call FFT(BxR(0,0),BxR(1,0),2,n,-1)
      Do 40 x=0,n-1
         Bx(x)=Bx(x)/(n*1.0d0)
 40   Continue

      Do 50 x=0,n-1
         Phix(x)=-(cdabs(B(x)))**2
 50   Continue

      Call FFT(PhixR(0,0),PhixR(1,0),2,n,1)
      Do 60 i=0,n/2
         Phix(i)=-0.50d0*i/s*Phix(i)
 60   Continue
      Do 70 i=n/2+1,n-1
         Phix(i)=-0.50d0*(n-i)/s*Phix(i)
 70   Continue
      Call FFT(PhixR(0,0),PhixR(1,0),2,n,-1)
      Do 80 x=0,n-1
         Phix(x)=Phix(x)/(n*1.0d0)
 80   Continue

      Do 90 i=0,n-1
          k(i)=mu2*ii*(cdabs(B(i)))**2*B(i)
          k(i)=k(i)-mu5*epsilon*(cdabs(B(i)))**2*Bx(i)
          k(i)=k(i)+mu6*2.0d0*epsilon*ii*B(i)*phix(i)
 90   Continue

      End



      Subroutine sixthorder(q,s,h,k0,epsilon,delta,gamma)
C***********************************************************************
C     Sixth-order operator splitting.
C***********************************************************************
      Complex*16 A
      Real*8 s,h,k0,epsilon,l2n,avewn,delta,qnum,w3,w2,w1,w0,hd2,AR
      Real*8 gamma
      Integer q,j,n
      Parameter(n=1024)
      Parameter(w3=0.7845136104775600d0)
      Parameter(w2=0.2355732133593570d0)
      Parameter(w1=-1.177679984178870d0)
      Parameter(w0=1.3151863206839060d0)
      Dimension A(0:n-1),AR(0:1,0:n-1)
      Equivalence(A,AR)
      Common A

      hd2=h/2.0d0

      Call FFT(AR(0,0),AR(1,0),2,n,1)
      write(2,11,advance="no")(cdabs(A(40-i))/(n*1.0d0),i=1,40)
      write(2,11) (cdabs(A(n-i))/(n*1.0d0),i=1,39)
 11   format(' ',79f9.6)
      Call FFT(AR(0,0),AR(1,0),2,n,-1)
      Do 30 i=0,n-1
         A(i)=A(i)/(1.0d0*n)
 30   Continue
      Call l2norm(l2n)
      Call awn(s,avewn)
      Call qnorm(s,qnum)

      Write(3,5) l2n,avewn,qnum
      Print *,'L2-norm, AWN, and Q',l2n,avewn,qnum
 5    Format(f22.16,f22.16,f22.16)

      Call linear(s,0.50d0*w3*hd2,k0,delta,epsilon,gamma)

      Do 10 j=1,q-1
         Call nonlinear(0.50d0*w3*h,k0,epsilon,s,gamma)
         Call linear(s,0.50d0*(w3+w2)*hd2,k0,delta,epsilon,gamma)
         Call nonlinear(0.50d0*w2*h,k0,epsilon,s,gamma)
         Call linear(s,0.50d0*(w2+w1)*hd2,k0,delta,epsilon,gamma)
         Call nonlinear(0.50d0*w1*h,k0,epsilon,s,gamma)
         Call linear(s,0.50d0*(w1+w0)*hd2,k0,delta,epsilon,gamma)
         Call nonlinear(0.50d0*w0*h,k0,epsilon,s,gamma)
         Call linear(s,0.50d0*(w0+w1)*hd2,k0,delta,epsilon,gamma)
         Call nonlinear(0.50d0*w1*h,k0,epsilon,s,gamma)
         Call linear(s,0.50d0*(w1+w2)*hd2,k0,delta,epsilon,gamma)
         Call nonlinear(0.50d0*w2*h,k0,epsilon,s,gamma)
         Call linear(s,0.50d0*(w2+w3)*hd2,k0,delta,epsilon,gamma)
         Call nonlinear(0.50d0*w3*h,k0,epsilon,s,gamma)
         Call linear(s,0.50d0*w3*h,k0,delta,epsilon,gamma)
         Call FFT(AR(0,0),AR(1,0),2,n,1)
         write(2,11,advance="no")(cdabs(A(40-i))/(n*1.0d0),i=1,40)
         write(2,11) (cdabs(A(n-i))/(n*1.0d0),i=1,39)
         Call FFT(AR(0,0),AR(1,0),2,n,-1)
         Do 40 i=0,n-1
            A(i)=A(i)/(1.0d0*n)
 40      Continue
         Call l2norm(l2n)
         Call awn(s,avewn)
         Call qnorm(s,qnum)
         Write(3,5) l2n,avewn,qnum
         
 10   Continue

      Call nonlinear(0.50d0*w3*h,k0,epsilon,s,gamma)
      Call linear(s,0.50d0*(w3+w2)*hd2,k0,delta,epsilon,gamma)
      Call nonlinear(0.50d0*w2*h,k0,epsilon,s,gamma)
      Call linear(s,0.50d0*(w2+w1)*hd2,k0,delta,epsilon,gamma)
      Call nonlinear(0.50d0*w1*h,k0,epsilon,s,gamma)
      Call linear(s,0.50d0*(w1+w0)*hd2,k0,delta,epsilon,gamma)
      Call nonlinear(0.50d0*w0*h,k0,epsilon,s,gamma)
      Call linear(s,0.50d0*(w0+w1)*hd2,k0,delta,epsilon,gamma)
      Call nonlinear(0.50d0*w1*h,k0,epsilon,s,gamma)
      Call linear(s,0.50d0*(w1+w2)*hd2,k0,delta,epsilon,gamma)
      Call nonlinear(0.50d0*w2*h,k0,epsilon,s,gamma)
      Call linear(s,0.50d0*(w2+w3)*hd2,k0,delta,epsilon,gamma)
      Call nonlinear(0.50d0*w3*h,k0,epsilon,s,gamma)
      Call linear(s,0.50d0*w3*hd2,k0,delta,epsilon,gamma)

      Call FFT(AR(0,0),AR(1,0),2,n,1)
      write(2,11,advance="no")(cdabs(A(40-i))/(n*1.0d0),i=1,40)
      write(2,11) (cdabs(A(n-i))/(n*1.0d0),i=1,39)
      Call FFT(AR(0,0),AR(1,0),2,n,-1)
      Do 60 i=0,n-1
         A(i)=A(i)/(1.0d0*n)
 60   Continue

      Call l2norm(l2n)
      Call awn(s,avewn)
      Call qnorm(s,qnum)
      Write(3,5) l2n,avewn,qnum
      Print *,'L2-norm, AWN, and Q',l2n,avewn,qnum

      End



      Subroutine l2norm(l2n)
C***********************************************************************
C     This subroutine calculates the M norm.
C***********************************************************************

      Complex*16 A
      Real*8 l2n
      Integer n,i
      Parameter (n=1024)
      Dimension A(0:n-1)
      Common A

      l2n=0.0d0
      Do 10 i=0,n-1
         l2n=l2n+(cdabs(A(i)))**2
 10   Continue

      l2n=l2n/(n*1.0d0)

      Return
      End

      Subroutine awn(s,avewn)
C***********************************************************************
C     This subroutine calculates the P norm.
C***********************************************************************
      Complex*16 A,B,ii,awnn
      Real*8 pi,s,BR,avewn
      Integer n,x,i
      Parameter(n=1024)
      Dimension A(0:n-1),B(0:n-1),BR(0:1,0:n-1)
      Common A
      Equivalence(B,BR)
      pi=4.0d0*datan(1.0d0)
      ii=dcmplx(0.0d0,1.0d0)

      Do 10 x=0,n-1
         B(x)=A(x)
 10   Continue

      Call FFT(BR(0,0),BR(1,0),2,n,1)
      
      Do 20 i=0,n/2
         B(i)=-ii*i/s*B(i)
 20   Continue
      Do 30 i=n/2+1,n-1
         B(i)=ii*(n-i)/s*B(i)
 30   Continue

      Call FFT(BR(0,0),BR(1,0),2,n,-1)

      Do 40 x=0,n-1
         B(x)=B(x)/(n*1.0d0)
 40   Continue

      awnn=0.0d0

      Do 50 x=0,n-1
         awnn=awnn+A(x)*dconjg(B(x))-B(x)*dconjg(A(x))
 50   Continue

      awnn=0.50d0*ii*awnn/(n*1.0d0)

      avewn=dreal(awnn)

      Return
      End

      Subroutine qnorm(s,qnum)
C***********************************************************************
C     This subroutine calculates the Q norm.
C***********************************************************************
      Complex*16 A,B,ii,qnumm
      Real*8 pi,s,BR,qnum
      Integer n,x,i
      Parameter(n=1024)
      Dimension A(0:n-1),B(0:n-1),BR(0:1,0:n-1)
      Common A
      Equivalence(B,BR)
      pi=4.0d0*datan(1.0d0)
      ii=dcmplx(0.0d0,1.0d0)

      Do 10 x=0,n-1
         B(x)=A(x)
 10   Continue

      Call FFT(BR(0,0),BR(1,0),2,n,1)
      
      Do 20 i=0,n/2
         B(i)=-ii*i/s*B(i)
 20   Continue
      Do 30 i=n/2+1,n-1
         B(i)=ii*(n-i)/s*B(i)
 30   Continue

      Call FFT(BR(0,0),BR(1,0),2,n,-1)

      Do 40 x=0,n-1
         B(x)=B(x)/(n*1.0d0)
 40   Continue

      qnumm=0.0d0
      Do 50 x=0,n-1
         qnumm=qnumm+B(x)*dconjg(B(x))
 50   Continue
      qnum=dreal(qnumm/(n*1.0d0))

      Return
      End

C-----------------------------------------------------------------------
      SUBROUTINE FFT (A,B,IS,N,ID)
C-- +------------------------------------------------------------------+
C-- | A CALL TO  FFT  REPLACES THE COMPLEX DATA VALUES  A(J) + i B(J), |
C-- | J=0,1,..,N-1 WHITH THEIR TRANSFORM                               |
C-- |                                                                  |
C-- |                             2 i ID PI K J / N                    |
C-- |    SUM     (A(K) + i B(K)) e                   ,  J=0,1,..,N-1   |
C-- |  K=0..N-1                                                        |
C-- |                                                                  |
C-- | INPUT AND OUTPUT PARAMETERS                                      |
C-- |     A   ARRAY A(0:*), REAL PART OF DATA/TRANSFORM                |
C-- |     B   ARRAY B(0:*), IMAGINARY PART OF DATA/TRANSFORM           |
C-- | INPUT PARAMETERS                                                 |
C-- |    IS   SPACING IN MEMORY BETWEEN CONSECUTIVE ELEMNTS IN A AND B.|
C---|         (USE IS=+1 FOR ELEMENTS STORED CONSECUTIVELY)            |
C-- |     N   NUMBER OF DATA VALUES, MUST BE A POWER OF TWO            |
C-- |    ID   USE  +1 OR -1  FOR DIRECT AND INVERSE TRANSFORM RESP.    |
C-- |                                                                  |
C-- +------------------------------------------------------------------+
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION A(0:*),B(0:*)
      J = 0
C---   APPLY PERMUTATION MATRIX   --------
      DO 20 I=0,(N-2)*IS,IS
         IF (I.LT.J) THEN
            TR   = A(J)
            A(J) = A(I)
            A(I) = TR
            TI   = B(J)
            B(J) = B(I)
            B(I) = TI
         ENDIF
         K = IS*N/2
   10    IF (K.LE.J) THEN
            J = J-K
            K = K/2
            GOTO 10
         ENDIF
   20    J = J+K
C---   PERFORM THE  LOG2 N  MATRIX-VECTOR MULT.  ---
      S =  0.0D0
      C = -1.0D0
      L = IS
   30    LH = L
         L  = L+L
         UR = 1.0D0
         UI = 0.0D0
         DO 50 J=0,LH-IS,IS
            DO 40 I=J,(N-1)*IS,L
               IP = I+LH
               TR = A(IP)*UR-B(IP)*UI
               TI = A(IP)*UI+B(IP)*UR
               A(IP) = A(I)-TR
               B(IP) = B(I)-TI
               A(I) = A(I)+TR
   40          B(I) = B(I)+TI
            TI = UR*S+UI*C
            UR = UR*C-UI*S
   50       UI = TI
         S = DSQRT (0.5D0*(1.0D0-C))*ID
         C = DSQRT (0.5D0*(1.0D0+C))
         IF (L.LT.N*IS) GOTO 30
      RETURN
      END
